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SPECTRAL  REPRESENTATION  OF  TRANSIENT  WAVEFORMS 

by 

Charles  R.  Arnold 
Honeywell  Radiation  Center 


ABSTRACT 


For  steady-state  waveform  analysis,  the  classical  (possibly  smoothed) 
periodogram  of  the  sampled  waveform  gives  one  an  adequate  spectral 
representation.  For  transient  waveforms  of  unknown  duration  in  noise, 
however,  the  periodogram  generally  fails  in  that  it  is  tied  to  a fixed 
time  interval.  As  an  alternative,  a digital  computer  program  has  been 
developed  which  will  do  time-varying  spectral  estimation. 


Briefly,  the  program  may  be  described  as  a digital  equivalent  of  a con- 
stant Q comb  filter  bank  wherein  one  can  vary  the  frequency  range  covered 
and  the  frequency  resolution  (i.  e.  the  Q).  For  a given  specified  frequency 
range,  as  one  increases  the  frequency  resolution,  the  program  automatic- 
ally selects  more  filters  and  spaces  them  so  as  to  cover  the  specified 
frequency  range;  the  various  contiguous  filters  being  overlapped  at  the 
-3  dB  points.  The  instantaneous  energy  contained  in  each  filter  is  used  to 
modulate  the  z-axis  of  a CRT  display  and  hence  provide  a time-frequency- 
intensity  plot  of  the  time  varying  spectrum. 


I.  INTRODUCTION 

Consider  the  following  problem:  Given  sample  functions  of  various  tran- 
sient waveforms  embedded  in  noise,  how  may  one  determine  whether 
there  is  any  natural  grouping  of  the  transients  into  classes  causally 
related  to  the  mechanism  which  generated  the  transients.  From  past 
work  in  steady -state  system  evaluation,  one  knows  that  often  the  power 
spectrum  of  the  data  is  a much  more  consistent  basis  for  classification 
than  the  raw  noise  corrupted  waveform.  Thus,  one  is  motivated  to 
determine  a spectral  representation  (or  signature)  of  the  transients. 

As  to  practical  realization  of  the  above  problem,  one  may  think  in  terms 
of: 

1.  The  classification  of  passing  vehicular  traffic  in  a 
noisy  acoustic  background 

2.  The  classification  of  both  man-made  and  biologically 
generated  transients  in  an  ocean  background 

3.  The  discrimination  between  seismic  events  (earth- 
quake/blast) in  a noisy  background 


Although  motivated  by  these  classification  problems,  this  paper  is 
addressed  only  to  the  problem  of  obtaining  a spectra  representative  of 
the  transient  waveforms.  Additional  details  on  the  classification  problem 
may  be  found  in  the  accompanying  paper  by  Swanlund. 

If  the  transient  waveforms  presented  for  classification  are  well  above  the 

noise  to  the  extent  that  their  epoch  and  duration  are  easily  determined, 

then  any  of  the  methods  summarized  in  the  recent  paper  "Burst 

2 

Measurements  in  the  Frequency  Domain"  by  Cochran  et  al  , may  be  used 
to  generate  a spectral  pattern.  However,  the  transient  waveforms  are 
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usually  so  far  down  into  :he  noise  that  their  initial  epoch  and  duration  are 
not  readily  discernible.  Moreover,  one  is  often  faced  with  multipath 
problems  and  multiple  arrivals  of  the  transient. 

To  illustrate  these  problems,  consider  an  example  from  the  seismic  area. 
Figure  1 presents  a record  of  approximately  four  minutes  of  the  output  of 
a single  seismometer.  To  a trained  eye,  the  trace  does  contain  a seismic 
event  which  happens  to  be  a moderately  strong  earthquake.  The  initial 
epoch  of  this  seismic  transient  is  only  approximately  discernible  and  its 
duration  is  obscurred  by  multiple  arrivals.  There  are  many  other  seismic 
events  of  interest  which  are  even  much  further  down  into  the  noise  and  even 
their  occurrence  is  not  discernible  from  a single  trace. 


Figure  1 SEISMOGRAM  CONTAINING  EARTHQUAKE  EVENT 


2 

All  of  the  methods  of  Cochran  et  al,  as  presented  fail  when  the  transient’s 
epoch  and  duration  are  not  known.  With  the  duration  of  the  transient  un- 
known, one  may  conceivably  choose  some  upper  bound  Tmax  as  the  duration 
of  any  transient  of  interest  and  then  segment  one  s data  into  intervals  of 
length  Tmax  and  apply  any  one  of  the  techniques  of  reiei  ence  2 to  each  data 
segment. 
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This  is  not  satisfactory,  in  that  this  arbitrary  segmentation  ot  the  data 
may  cut  some  occurrences  in  half  and  thus  distort  and/or  lose  the  true 
spectral  pattern.  Next,  one  may  apply  the  methods  of  reference  2 to 
overlapping  time  intervals  of  the  data  but  this  leads  to  excessive  equip- 
ment requirements  in  the  case  of  analogue  signals  and  excessive  compu- 
tational requirements  in  the  case  of  sampled  data. 

As  an  alternative,  this  paper  presents  a technique  which  has  been  im- 
plemented upon  a digital  computer  and  which  provides  a time- varying 
spectral  estimate.  Briefly,  the  program  may  be  described  as  a digital 
equivalent  of  a constant  Q comb  filter  bank  wherein  one  can  vary  the  fre- 
quency range  covered  and  the  frequency  resolution  (i.  e. , the  Q) . This 
program  has  been  found  to  be  most  effective  against  transients  in  that 
it  is  not  tied  to  a fixed  interval  or  initial  epoch.  Moreover  , from  re- 
sults to  be  developed  below,  the  program  does  not  require  the  detectors 
and  integrators  (low-pass  filters)  usually  associated  with  filter  banks  (cf. 
Figure  2 of  Reference  2)  and  hence  does  not  require  the  additional  response 
time  which  degrades  transient  analysis.  A complete  description  of  the 
program  is,  however,  postponed  until  after  the  algorithm  it  uses  has  been 
motivated  and  mathematically  justified.  This  is  not  because  of  the  com- 
plexity of  the  algorithm  but  just  the  opposite,  in  that  its  simplicity  is  best 
appreciated  only  after  having  gone  through  its  motivation. 

Thus,  in  the  next  section  begins  a sequence  of  definitions,  results  and 
related  discussions  most  of  which  are  well  known  and  well  worn.  Some 
are  hopefully  new,  but  if  not  their  collection  here  may  be  justified  in 
that  they  represent  the  path  down  which  this  author's  thinking  has  evolved. 
The  end  result  being  the  computer  program  which  has  been  found  to  be  a 
very  eftective  tool  in  many  areas  of  application. 
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II.  THE  PERIODOGRAM 


Given  a sample  of  continuous  time  data  X(t)  of  duration  T,  the  classical 
Schuster  periodogram  (or  sample  spectral  density  function)  is  defined  by 


1-p  (u>)  = -TjT 


X(t)  ej"ldt 


(2.1) 


Similarly,  for  sampled  data  where  the  sample  consists  of  N values 
equally  spaced  in  time,  one  uses 

N 


*N  <">  * T 


X e 
n 


jcun 


n=l 


(22) 


In  either  case,  if  one  chooses  a discrete  set  of  frequencies,  say 


wk  = 


27rk 

~T~ 


or  «k  = 


277k 

IT 


(2.3) 


and  expands  the  complex  exponential,  one  finds,  for  example,  that 

2 


!T  (wk)  = T { 


/ 

T 

r 

j X(t)cosjl^ 

dt 

✓ 

0 

(2.4) 


X (t)  sin  dt 
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*Wr* 


from  which  one  can  recognize  that 


IT(wk)  = T(a£  + b£)  (2.5) 

where  the  a^  and  are  the  usual  Fourier  coefficients.  Thus  it  is  quite 
natural  to  expect  that  when  the  data  consists  of  a truly  periodic  signal 
well  above  the  noise,  as  is  often  the  case  for  steady-state  phenomena, 
the  periodogram  provides  a completely  effective  analytical  tool. 


When  one  is  interested  in  X (t)  as  a stochastic  process,  one  naturally  turns 
to  consideration  of  the  sample  auto-covariance  function  as  a vehicle  for 
spectral  estimation  since  theoretically  the  spectral  density  P(ai)  is  nothing 
more  than  the  cosine  transform  of  the  auto-correlation  function  of  the  pro- 
cess. Specifically  for  discrete-time  data,  one  introduces  the  sample  auto- 
covariance sequence. 


N-v 


n=l 


X 


n+v 


from  which  it  follows  that 

W = Co  + 2 


N-l 

V” 


n=l 


Cn  cos  (nwk) 


(2.6) 


(2.7) 


The  validity  of  (2.  7)  is  a direct  consequence  of  algebraic  identities  and  is 
independent  of  any  assumptions  about  the  nature  of  the  process  X,  This 
fact  is  often  lost  by  some  people  to  the  extent  that  they  will  produce  two 
different  computer  programs.  The  first  program  employing  the  direct 
approach  (2  2)  to  be  used  on  deterministic  data  and  a second  program  using 
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(2.  6)  and  (2.  7)  to  be  used  upon  random  or  noisy  data.  The  truth  of  the 
matter  is,  given  N data  samples  X^,  both  programs  (at  best)  can  give 
N independent  values  of  the  sample  spectral  density  IN  (cu^)  and  these 
estimated  values  will  be  the  same  outside  of  minor  variations  introduced 
by  different  rounding  errors  in  the  two  programs. 

However,  as  is  well  known,  when  one  desires  a random  process  as  the 
underlying  model,  the  periodogram  has  its  problems.  It  is  true  that 
the  periodogram  is  an  asymptotically  unbiased  estimate  of  the  true 
spectral  density  function  P(co),  i.  e. , 

Lim  E < I (o>)  \ ^ P (oj)  . (2. 8) 

T — 00  ( / 

The  difficulty  with  the  periodogram  is  that  when  used  against  a random 
process,  it  is  not  a consistent  estimator  in  that 

Lim  Var  |lT  (u>)  j — ^ P2  (u>).  (2.  9) 

This  lack  of  consistency  is  reflected  in  the  high  degree  of  variability 
between  the  spectral  estimates  from  successive  samples  of  the  same 
noise  process. 
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III. 


THE  SMOOTHED  PERIODOGRAM 
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In  1946,  Daniell  suggested  how  the  periodogram  could  be  made  into  a 
consistent  estimator  of  the  spectral  density.  His  solution  was  in  effect 
to  average  the  periodogram  over  neighboring  frequencies.  Since  then, 
there  has  evolved  the  general  concept  of  smoothing  the  periodogram  with 
some  spectral  window  W^,  (w)- 


The  smoothed  periodogram  is  given  by  the  convolution, 


PT  (w)  = 


WT  (w-A)  IT  (A)  dA, 


of  the  raw  periodogram  with  the  spectral  window. 


(3.1) 


In  the  last  20  years,  a multitude  of  smoothing  windows  have  made  their 
appearance.  Some  of  the  more  popular  ones  are  associated  with  the  names 
Bartlett,  Hamming,  Hanning,  Parzen,  etc.  There  is  still  a good  deal  of 
controversy  over  the  choice  of  any  specific  window  to  insure  an  optimal  es- 
timate, let  alone  any  initial  agreement  over  a criterion  of  optimality.  This 
paper  does  not  hope  to  explore  this  facet  of  the  problem  (the  interested  rea- 
der may  consult  references  4-9)  but  only  to  make  the  point  that,  generally, 
some  smoothing  is  required  when  the  periodogram  is  used  against  a random 
process. 


The  mechanization  of  (3. 1)  for  discrete-time  data  is  often  quite  simple  and 
clearly  illustrates  Daniell' s original  suggestion  of  averaging  the  periodogram 
over  neighboring  values.  A particular  case  being  the  discrete-time  Hanning 
smoothed  periodogram  which  is  given  by 
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Over  the  years,  this  author  has  found  this  smoothing  algorithm  to  be  as  ef- 
fective as  any  for  run-of-the-mill  problems.  Its  acceptance  by  this  author 
being  its  simplicity,  especially  for  implementation  on  a binary  computer. 


The  smoothed  periodogram  is  also  a reasonable  choice  as  an  estimator  of 
the  spectral  density  function  of  noise-free  transient  waveforms.  This  is 
because  the  smoothing  process  tends  to  remove  round-off  and  discretiza- 
tion errors  introduced  by  the  computation  and  yet  does  not  degrade  the  es- 
timate since  the  theoretical  spectrum  is  continuous  and  not  subject  to  sharp 

jumps  or  spikes.  As  an  example,  Figure  2 presents  the  smoothed  periodo- 

2 

gram  estimate  of  the  transfer  gain  characteristic  H (ju)  derived  from 
the  sampled  impulse  response  of  the  system, 

2 

H(s)  = - g (3'3) 

s + 2£  u)  s + w 
n n 

where 

d)n  = 6tt,  £ = 0.05,  At  = 0.005.  (3.4) 

Figure  3 presents  the  Hanning  smoothed  periodogram  for  the  complete  four- 
minute  seismogram  of  Figure  1.  The  earthquake  event  is  indicated  by  the 
double  spike  of  energy  about  1. 0 Hz  and  by  the  single  spike  at  2.  0 Hz.  The 
major  portion  of  the  energy  centered  about  0.2  Hz  is  due  to  the  micro- seis- 
mic noise.  Additional  discussion  is  postponed  until  Section  IX.  Specific 
details  of  the  program,  which  were  used  to  determine  the  estimates  of  Fig- 
ures 2 and  3,  may  be  found  in  reference  10. 
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IV. 


THE  FORCED  HARMONIC  OSCILLATOR 


L 


A.  The  Continuous-Time  Case 

Consider  the  forced  harmonic  oscillator  given  in  state- vector  form  by 

V 


xk(t)  = Ak  xR  (t)  + b u(t)  = 


-“k  0 


xlk(ty  ' o 

lk  ' + , u(t). 

I1] 


-x2k(t) 


(4.1) 


The  fundamental  matrix  for  this  system  is  readily  found  to  be 
A,_t  cos  (w,  t)  -sin(u).t)' 


\(t)  = e 


sin  (cu'^t)  cos  (aJkt)_ 


(4.2) 


and  the  general  solution  of  (4. 1)  for  xk(0)  = 0 is  given  by 


11 


r a (T-t)  ~k 

x^(T)=  e b u(t)dt  = e 


a,  t r 


, cos  (w,  t)  sin  (w.t)  0 

I k k u(t)dt  (4.3) 

^0  -sin  (wkt)  cos  (wkt)_  _1_ 


Now  consider  the  norm  of  the  system  (4. 1)  or  more  specifically,  the  quantity 


_1 

T 


*k  (T) 


AkT 


By  inspection  of  (4.2),  one  can  see  that  e is  an  orthogonal  matrix  and 

AkT 

is  thus  norm-preserving  (i.  e. , e is  a pure  rotation).  Hence,  from  (4.3) 

i:2 


one  has 
1 


X(T) 


_1 

T 


«f 


u(t)  sin  (o.'kt)  dt 


I!  T 


j ^ u(t)  cos  (wkt)  dt 


^ jjXjjT)  !;  = u(t)  cos  (wRt)  dt  +/  u(t)  sin  (u>kt)  dt  (4.4) 

" J ‘ rO  KO  I 


A comparison  of  (4.  4)  with  (2.4)  now  yields  a principal  result  for  this  paper, 
namely, 


T 


- IT  (wk) 


(4.5) 


which  shows  that  the  classical  periodogram  for  a discrete  set  of  frequencies 

can  be  determined  by  driving  a set  of  undamped  harmonic  oscillators  (initially 

1 

quiescent)  with  the  input  data  sample  and  then  measuring  ^ times  the  square 
of  the  norms  of  the  oscillator's  state- vector. 


B.  The  Discrete-Time  Case 

Rather  than  start  from  a difference  equation,  this  paper  will  first  derive  it 
from  the  continuous  time  case  in  order  to  motivate  the  results  of  this  section 
plus  other  results  below.  Thus,  returning  to  an  equivalent  form  of  (4. 3)  for 
the  continuous- time  case,  one  can  write 


fT  V 

xk(T)  = e b u (T-t)dt. 


(4.6) 


Now  introduce  the  usual  simple  minded  discrete-time  approximation 
(t  — nAt,  T -NAt,  u(t)  -u(n),  etc.), 


I 


rrk  (N)  = At  ) B“  bu(N-n) 
n=0 


where 


Bk  * e 


V‘ 


From  (4.  7a),  one  can  write 

N 


B,n  b u(N  + 1 - n) 


xfe(N  + 1)  = At  x^k  « 

n=0 
N 


= At  B^  b u(N  + 1 - n)  + At  lb  u(N  + 1) 

“n=l 

N-l 

= At  B^+1  b u(N  - X)  + At  b u(N  + 1) 

x=o 


N-l 

= Bk  At  ^ Bk  b u(N  - X)  + At  b u(N  + 1) 
X=0 


xk  (N  + 1)  = Bk  sk(N)  + At  b u(N  + 1) 

Equation  (4.8a)  which  can  be  written  in  greater  detail  as 
xik  (N  + !)  | | cos  (u>kAt)  -sin  (a>kAt) 


x2k  + 1 ) 


xlk(N)' 

V 

+ 

■J 

*2k<N) 

1 

At  u(N  + 1) 


yields  an  algorithm  for  discrete  filters  of  a comb  filter  program. 


(4.7a) 


(4.7b) 


(4. 8a) 


(4.8b) 
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n 


Now  returning  to  the  equivalence  of  the  discrete- time  forced  harmonic 

oscillator  (4.8)  and  the  discrete-time  periodogram,  consider  again  for 
277k  _ 27Tk 
NAt  ’ 


wk  T 


Bk  = e 


A.  At 
k 


/27Tkv 
cos  ( N ) 

, 277k  ■ 
-sin  (^-, 

. /277k^ 

sm  (-jj-) 

,277k 
cos  (^- 

By  induction,  one  can  readily  show  that 


, 27rkn  . 
cos  ( N ) 

. , 2tfkn  v 
-sm  ( N ) 

Bk  ■ 

1 n 

. / 2i7kn  x 

Lsin(  N ) 

/ 277knv 
cos  ( N ) 

1 * * 

From  (4.  9)  and  (4. 10),  it  follows  that 


= Bk’  Bk  = h B?"n  = B^  = B" 


and 


Rn 

Bk 


, 277 kn  \ 
cos  (-Jj-) 

/ 277kn 

sm  ( N 

, 277knx 
[_-sin(  N ) 

/ 277kn 
cos  ( N 

(4.9) 


(4.10) 


(4.11) 


(4.12) 


n 


V. 


DISCUSSION 


At  this  point  let  us  re-examine  the  total  problem  of  spectral  estimation. 
This  may  be  summarized  by  the  following  chart: 


Process 

Spectrum 

Estimator 

Periodic 

Stochastic 

Transient 

Discrete 

Continuous 

Continuous 

Periodogram 

Smoothed  Periodogram 
Smoothed  Periodogram 

For  most  real-world  problems,  however,  one  does  not  have  isolated  proces- 
ses, but  usually  a combination,  such  as  periodic  and  transient  waveforms  in 
a noise  background.  Yet,  one  is  interested  in  determining  a single  estimator 
to  work  against  the  combined  process.  From  what  has  been  offered  to  date, 
one  is  reasonably  justified  in  using  a smoothed  periodogram  in  all  cases,  and 
realizing  that  discrete  spectral  lines  are  spread  somewhat  by  the  estimator. 
The  smoothed  periodogram  is  a distinct  possibility  for  the  case  of  interest 
in  this  paper,  namely  for  transient  in  a noise  background.  So  far,  this  paper 
has  offered  three  computation  techniques  for  determining  the  periodogram 
which  may  be  summarized  briefly  as: 


T-l.  The  direct  method  via  (2. 2). 

T-2.  The  correlation  method  via  (2. 6)  and  (2. 7). 

T-3.  The  bank  of  forced  harmonic  oscillators  given  by 
(4. 8)  with  (4. 15). 


f 
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Some  notable  difficulties  with  the  periodogram  for  the  problem  at  hand  are, 
however,  summarized  as  follows: 

Dl.  The  periodogram  must  be  smoothed  in  order  to  obtain  a con- 
tinuous spectra  for  transients  and  a consistent  estimator  against 
noise. 

D2.  The  periodogram  is  tied  to  a fixed  time  interval  hence,  both  the 
epoch  and  duration  of  any  transient  must  be  known. 

D3.  The  periodogram  requires  overlapping  estimate  to  be  made  upon 
evolutionary  processes  in  order  to  determine  the  time- varying 
spectrum. 

The  next  section  presents  a variation  on  technique  T-3  above,  which  meets 
all  the  stated  objectives  to  the  periodogram  and,  moreover,  the  new  techni- 
que also  circumvents  a major  difficulty  of  other  techniques  as  offered  by 
2 

Cochran  et  al,  namely, 

D4.  Other  techniques  require  detection  and  integration  with  subsequent 
loss  in  response  time  crucial  to  transient  analysis. 
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VI. 


THE  FORCED  DAMPED  HARMONIC  OSCILLATOR 


This  section  presents  the  modifications  resulting  from  the  introduction  of 
some  damping  into  the  harmonic  oscillators  of  Section  IV  and  then  goes  on 
to  consider  how  this  modification  to  the  technique  T-3  yields  a very  effi- 
cient algorithm  for  spectral  estimation  free  of  all  the  objections  raised  in 
the  previous  section. 

In  order  to  add  some  damping  to  the  formulation  for  the  forced  harmonic 
oscillator  of  Section  IV,  one  can  make  the  following  change  in  th^  transi- 
tion rririx  of  the  system. 


0 

"wk 

r 0 

-s 

-wk 

0 

^ k 

.wk 

-«V 

Another  choice  and  the  one  used  by  the  author  is: 


"-a  - <jj  ' 

A = 1 ° 

Ak 


1% 


- a 


r-H 

-H 


For  this  particular  choice,  the  new  fundamental  matrix  for  the  k-th  filter 
is  readily  found  to  be 


r -at  / *\  -ai 

A.  t e cos  K'1  _e 

* = e k = 

-at  , / .v  -at 
e sin  (^0w  e 

where  as  usual  (£  is  called  the  damping  ratio) 
a = Cu>k  and  wQ  = «k"\/l  - £2 


sin  (u>  t) 
o 

cos  (u>  t) 
o J 


(6.1) 


(6.2) 


(6.3) 


(6.4) 
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The  discrete-time  algorithm  is  the  same  as  (4. 8a), 

\ (n+1)  = Bk  (n)  + At  b u(n+l),  (6.5) 

except  now  the  matrix  Bk  is  derived  from  (6. 3)  with  again 
AkAt 

Bk  = e (6.6) 

As  an  estimator  of  the  spectral  density  function,  this  author  proposes  that 
one  employ  a bank  of  filters  (damped  oscillators)  covering  the  frequency 
range  of  interest  and  using  the  algorithm  (6. 5).  By  analogy  with  the  un- 
damped case,  the  square  of  the  norms  of  the  recursion  algorithm's  state 
vector  is  again  taken  as  a spectral  measure  (see  Figure  4). 


Now  consider  how  this  proposed  estimator  meets  the  objections  raised  in 
the  previous  section.  First,  as  is  well  known,  the  introduction  of  damp- 
ing detunes  the  oscillator.  This  detuning  of  the  harmonic  oscillator  has 
the  effect  that  it  widens  the  bandwidth  of  the  system  when  considered  as 
a bandpass  filter.  But  this  is  just  the  desired  effect  for  the  smoothing  of 
the  periodogram  (Dl). 

Secondly,  the  filter  algorithms  can  be  run  ''open-loop, " since  the  damping 
automatically  throws  away  the  oldest  past  input.  Moreover,  because  of  its 
recursive  nature,  the  algorithm  is  quite  efficient  numerically.  This  non- 
reset capability,  coupled  with  its  high  numerical  efficiency,  allows  one  to 
start  the  algorithms  well  ahead  of  the  initial  epoch  of  any  transient.  The 
filter  may  then  be  stepped  for  any  desired  duration  and  numerous  over- 
lapping spectral  estimates  may  be  obtained  at  will  by  subsequently  calcula- 
ting the  norms  of  the  algorithms  state  vectors.  Hence  the  second  and  third 
objection  (D2  and  D3)  are  readily  met.  Finally,  the  use  of  the  norms  as  a 
spectral  measure  circumvents  the  loss  of  response  time  due  to  the  usual 
detection  and  integration  of  other  techniques  (D4). 
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ADDITIONAL  ALGORITHMS  FOR  THE  COMPUTER  PROGRAM 


A consideration  of  the  input/output  transfer  functions  for  the  two  components 
of  the  damped  harmonic  oscillator’s  state- vector  shows  them  to  be  a linear 
combination  of  the  following: 

2 


Hl(s)  = -j 


s + 2£o>,  s + uv 


Hgts)  = 


S GO, 


s + 2Co>k  s + o>k 


(7.1) 


From  Figure  2,  one  can  see  that  the  first  of  these  H^(s)  has  significant 
gain  at  the  low  end  of  the  spectrum.  This  response  to  d-c  and  other  low 
frequencies  by  the  periodogram  is  well  known  and  is  usually  circumvented 
in  other  techniques  by  various  mean-removal  and  detrending  operations  prior 
to  computation  of  the  periodogram.  As  an  alternative  to  a detrending  opera- 
tion, the  present  program  differentiates  the  input  time  series  before  applica- 
tion to  the  various  filter  algorithms.  For  discrete-time  data,  this  is  nothing 
more  than  using  the  first  differences 


u(n+l)  = Xn+1  - Xn 


(7.2) 


to  drive  the  filter  algorithms  (6.  5). 


With  the  introduction  of  damping  into  the  second  order  algorithms,  one  is 

12 

restricted  to  a finite  bandwidth  Af  given  by  the  usual  relation 


q = Ji L = f 

w Aw 


Af 


1 

2? 


(7.3) 


Also,  for  many  applications,  a log  frequency  scale  for  the  estimated  spectrum 
is  desired.  This  subsequently  dictated  that  the  bank  of  filters  be  derived  on  a 
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basis  of  constant  quality  Q.  The  usual  procedure  is  to  overlap  the  various 
contiguous  filters  at  their  -3dB  points  as  in  Figure  5.  For  any  given  Q or 
equivalently  for  any  given  damping  ratio  C,  the  discrete  set  of  center  fre- 
quencies for  the  filter  bank  are  readily  derived  from  the  algorithm 


J 
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Figure  5 TRANSFER  GAIN  CHARACTERISTICS 
OF  COMB  FILTER  BANK 
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For  a damping  ratio  of  £ = 0. 115  and  a base  frequency  of  2.  5 Hz,  this  algo- 
rithm produces  a program  which  makes  the  computer  look  like  a standard 
1/3-octave  analyzer  with  the  set  of  center  frequencies  as  given  in  Table  1. 


Table  1 


2.79 

28.1 

283 

3.51 

35.4 

357 

4.  43 

44.6 

449 

5.58 

56.2 

566 

7.02 

70.8 

713 

8.85 

89.2 

898 

11.15 

112.3 

1132 

14.05 

141.5 

1426 

17.70 

178.3 

1797 

22.30 

224.7 

2264 

Although  the  proposed  estimation  scheme  does  not  require  any  post-detection 
smoothing,  some  time  smoothing  of  the  state-vector  norms  has  been  found 
beneficial,  especially  for  display  purposes.  To  date,  two  modes  of  smooth- 
ing have  been  implemented.  The  first  smoothes  (in  time)  each  output  vector 
norm  with  a low-pass  filter  whose  time  constant  is  fixed,  and  the  same  for 
all  frequency  bands.  The  second  mode  smoothes  each  output  norm  by  a low- 
pass  filter  whose  time  constant  is  some  fixed  proportion  of  the  respective 
I filter's  time  constant.  In  either  the  fixed  or  proportionally  smoothed  case, 

the  low- pass  filter  algorithm  is  given  by 


Pk(n+1)  = bkPk(n)+ak  xk(n+l)  2,  (7.5) 


where  in  the  fixed  case, 

ak  = Q0At,  bk  = exp(-aAt);  (7.6a) 

and  in  the  proportional  case 

ak  = ^27rfkNAt  At’  bk  = exP  ("ak)  (7.6b) 

for  some  choice  of  the  parameter  /3.  (The  fixed  point  algorithm  of  the  pro- 
gram requires  that  /3N^t  = 0.  3. ) 


VII!. 


THE  COMPUTER  PROGRAM 


The  algorithms  proposed  in  this  paper  have  been  implemented  upon  a 
Honeywell,  Computer  Control  Division's  DDP-24  general  purpose  com- 
puter. The  skeleton  of  the  computer  program  is  coded  in  FORTRAN, 
but  all  the  plotting  and  repetitive  calculations  are  coded  as  machine 
language  subroutines  for  greater  speed  and  efficiency.  Figure  6 pre- 
sents a flow  chart  of  the  complete  program.  The  program  is  presently 
configured  to  accept  the  input  data  as  records  of  500  samples  each  from 
magnetic  tape.  The  input  control  parameters  are: 

1.  NRCD:  The  number  of  input  data  records  to  be  processed. 

2.  NDT:  The  plot  (and  print)  increment. 

3.  DT:  The  basic  time  increment  At  of  the  data. 

4.  ZETA:  The  damping  ratio  f . 

5.  BETA:  The  output  norm  smoothing  ratio  /3  . 
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The  present  program  is  also  configured  to  cover  a three  decade  scale  in 
frequency  from  to  the  Nyquist  (or  folding)  frequency  fNy.  Since 


Ny 


1 

3~aT 


(8.1) 


the  program  can  thus  determine  via  algorithm  (7.  4)  with 


(8.2) 


'base  - <°-001>fNy 

the  number  of  the  filters  and  their  proper  spacing  so  as  to  cover  the  specified 
frequency  range.  With  a damping  ratio  of  £ = 0. 115  (Q  » 4.  3),  one  obtains 
thirty  1/3-octave  filter  as  denoted  in  Section  VIL  For  f = 0.03846,  ninety 
1 /9-octave  filters  are  obtained. 
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Figure  6 FLOW  CHART  OF  THE  COMB  FILTER  PROGRAM 
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Except  for  the  mode  of  output  display,  the  program  is  quite  standard. 
Fortunately,  the  Radiation  Center  has  a CRT  display  with  a 512  x 512 
x-y  resolution  interfaced  with  the  DDP-24.  In  addition,  the  z-axis  or 
intensity  of  the  display  has  16  levels  available  under  program  control. 

With  this  capability,  the  computer  program  is  configured  to  display  the 
evolution  of  a time -varying  spectrum.  Starting  from  the  top  of  the  dis- 
play, each  line  presents  an  instantaneous  measure  of  the  energy  contained 
in  each  of  the  algorithms  comprising  the  filter  bank,  (increasing  frequency 
left  to  right).  As  time  evolves  (and  more  data  is  processed)  the  horizontal 
trace  is  moved  down  the  face  of  the  CRT  to  provide  a time  axis  to  the 
energy  density  spectrum  (see  the  figures  in  the  next  section).  The  time 
scale  of  the  display  is  controlled  by  the  parameter  NDT  which  is  the 
multiple  of  the  basic  sample  spacing  At  (DT)  used  as  a display  increment; 
500  lines  comprising  one  frame.  The  intensity  of  the  display  for  any 
given  filter  (frequency  band)  is  actually  the  logarithm  (base  2)  of  the 
filter's  energy  and  hence  each  grey  scale  of  the  display  corresponds  to 
a 3 dB  difference  in  energy. 


The  entire  evolution  of  the  display  (and  thus  the  spectrum)  is  recorded  by 
a camera  whose  shutter  is  left  open.  After  one  frame  has  been  generated, 
the  computer  pauses  to  allow  the  operator  to  remove  and  reload  the  camera 
with  a new  Polaroid  film  pack.  Restarting  the  computer  provides  for  con- 
tinued evolution  of  the  data's  spectrum. 


IX.  SOME  RESULTS 


The  merits  of  the  comb  filter  technique  and  the  resulting  computer  pro- 
gram are  amply  demonstrated  by  the  results  presented  in  this  section. 

The  only  input  data  sample  considered  is  the  seismogram  of  Figure  1. 

The  data  sample  consisted  of  sampled  values  (at  20  samples  per  second) 
of  approximately  four  minutes  of  the  output  of  a single  seismometer.  The 
estimated  1 /9-octave  time  varying  spectrum  for  the  entire  data  sample  by 
the  comb  filter  program  is  presented  on  the  left  side  of  Figure  7.  On  the 
right  side  of  Figure  7 is  an  expanded  (in  time)  version  of  the  1/9  octave 
spectrum  of  the  central  portion  of  the  seismogram  which  contains  an 
earthquake.  Two  definite  arrivals  of  the  event  which  spectrally  is  a 
double  line  about  1.  0 Hz  plus  some  2.  0 Hz  energy  can  clearly  be  seen. 

Also,  the  after  effect  has  considerable  spectral  detail.  All  of  the  time- 
varying  detail  is  lost  by  a conventional  periodogram  analysis.  In  fact,  if 
all  the  energy  contained  in  each  frequency  band  could  be  projected  into  one 
time  plane,  one  would  recover  the  periodogram  estimate  of  the  total  sample 
as  given  in  Figure  3. 

Figure  8 duplicates  the  signature  on  the  right  in  Figure  7 plus  presenting 
an  enhanced  version  obtained  by  simply  restricting  the  CRT  intensity  to  a 
cruder  quantification.  In  both  Figure  7 and  8,  no  post-detection  (post- 
norm) smoothing  has  been  applied.  It  is  this  minimal  response  time  of 
the  computer  program  which  allows  all  the  time-varying  spectral  detail  to 
be  seen. 

Figures  9 and  10  present  signatures  comparable  to  those  of  Figure  7, 
except  in  Figure  9,  a uniform  (across  frequency)  post-detection  smooth- 
ing has  been  applied,  and  in  Figure  10,  the  post-detection  smoothing 
has  been  proportional  with  frequency.  In  both  cases  the  post -detection 
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Figure  9 


smoothing  tends  to  degrade  the  fine  time  structure  of  the  spectrum. 
However,  for  some  other  applications  of  the  program  wherein  the 
transients  had  narrow  band  signatures  of  longer  duration,  the  post- 
detection smoothing  helped  to  mitigate  the  effect  of  the  background 


X.  SUMMARY  AND  CONCLUSIONS 


This  paper  has  offered  an  alternative  technique  for  the  computation  of  the 
classical  periodogram  as  the  square  of  the  norms  of  the  state-vectors  of 
a set  of  undamped  harmonic  oscillators  forced  by  the  input  data  sample. 
Then  after  consideration  of  some  of  the  problems  with  the  periodogram, 
the  paper  has  shown  how  the  introduction  of  some  damping  into  the  pre- 
vious formulation  produces  extremely  efficient  algorithms  for  the  deter- 
mination of  time- varying  spectra. 

These  algorithms  have  been  implemented  into  a digital  computer  program 
which  has  subsequently  been  used  to  determine  the  spectral  signature  of 
transient  waveforms  embedded  in  noise  from  many  areas  of  application. 
The  high  numerical  efficiency  of  the  computer  program,  plus  the  high 
degree  of  detail  in  the  resulting  spectra  implies  that  the  technique  pre- 
sented in  this  paper  should  be  considered  in  any  future  spectral  estima- 
tion program. 
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